Inflow with Ram Pressure Stripping
The main version of the Inflow with Ram Pressure Stripping model (IRPS) and its gas-$\alpha$-abundance equations are here. I showed models for the range of [$\alpha$/H]$_{s}$, $\xi$ and $\zeta$ as examples.
This notebook is a simplified tutorial, from which you will learn:
-- how α-abundance distribution function (ADF) of IRPS model looks like
-- how this ADF uses gas as abundance-dependent function and how they are constructed
-- you will learn how to build this GCE-model for different parameters
numpy - used for numpy.arrays and conditional usage of indexes
scipy.special - the gas-metallicity equations has hypergeometrical and exponential functions in them
matplotlib.pyplot - for visualization
plotly - is required only for the last example with 3D figure (feel free to comment it)
In [1]:
import numpy as np
from scipy.special import hyp2f1, expi
import matplotlib.pyplot as plt
# for 3D figure (the very last example) only,
# You can comment these next two lines:
import plotly
import plotly.graph_objs as go
def format_ticks():
'''
This function is for figure formatting
'''
plt.rcParams.update({'font.family': 'DeJavu Serif'})
plt.rcParams.update({'font.serif':'Times New Roman'})
plt.gca().tick_params(which='major', direction="in", length=5.,\
labelbottom=True, labeltop=False, labelleft=True, labelright=False,
bottom=True, top=True, left=True, right=True)
plt.gca().tick_params(which='minor', direction="in", length=3.,\
labelbottom=True, labeltop=False, labelleft=True, labelright=False,
bottom=True, top=True, left=True, right=True)
Main functions
Part 1: gas--$\alpha$-abundance relations
This cell contains functions to build gas-abundance relation: eq_before_g, eq_after_g_abs, eq_after_g_equal build Z(g) before the forced gas loss occured (Z1(g) corresponds to "Before RPS" (ram pressure stripping)), and after (for ξ != ζ and ξ == ζ respectively; note, that ξ parameter I write everywhere as F, ζ parameter I write as E or Es; Z2(g) corresponds to "After RPS").
To merge them into a single solution, the initial conditions for equation(s) "After RPS" are replaced with RPS-value obtained from parameter [α/H]str and corresponding gas fraction point obtained with function g_rps.
long_coef is a part of the expression for Z2(g) ("After RPS") with ξ != ζ.
gas_equation combines above formulae and returns gas for equally spaced N = points alpha abundances within [α/H]min=-5, [α/H]max=0
There are solutions Z(g) which on g-log(Z) plot are convex. It depends on g-ξ, g-(ξ-ζ) signs: if gas is below than the normalized inflow, the function is monotonic, but if less, the ([α/H], g) plot will contain two "branches", upper, with gas decrease for α-abundance increase, and lower branch, where the gas decreases for α-abundance decrease. In this situation I chose "upper" branch, while fixation of inflow after [α/H]str for some minimal value also decreases the significance of the lower branch. Figure 1, green dashed line illustrates this point: the turning point is almost at g = 0.
In [2]:
def eq_before_g(F, z_inflow, z_initial, g, g_in = 1.):
'''
Determines z' = Z/p_eff = (10^[alpha/H])/p_eff
for array of g using equation "Before RPS"
(where p_eff = p/(1+eta) - effective yield)
------
Input:
------
F normalized inflow, xi
z_inflow alpha-abundance of inflow, expressed as z'_inflow
z_initial initial alpha-abundance, for z'_initial
g array of g values
g_in initial fraction of gas. We assumed that galaxy has started as cloud of gas: g_in=1
------
Returns:
------
z array of z' values
'''
z = np.full(len(g), -np.inf)
if F!=g_in:
i = np.where((((F-g)/(F-g_in)>0) & (g!=0.)))[0]
z[i] = g_in*(F-g[i])/g[i]/(F-g_in)*z_initial+\
(F-g[i])/g[i]*(np.log(((F-g[i])/(F-g_in)))) +\
F*(g[i]-g_in)/g[i]/(F-g_in)*(z_inflow+1)
return z
def g_rps(F, z_inflow, z_initial, z_str):
'''
Determines gas fraction in units of initial galaxies' mass
at moment when RPS commences, associated with [alpha/H]=[alpha/H]_s.
------
Input:
------
F normalized inflow, xi
z_inflow alpha-abundance of inflow, expressed as z'_inflow
z_initial initial alpha-abundance, for z'_initial
z_str alpha-abundance by achieving which RPS begins, z'_str
------
Returns:
------
g_str gas fraction at z' = z'_str = (10^[alpha/H]_s)/p_eff
'''
g = np.linspace(0., F*1.9999, 1000)
z = eq_before_g(F, z_inflow, z_initial, g)
z[np.where(np.isnan(z))] = -np.inf
delta = np.min(abs(z-z_str))
if z.max()>z_str:
g_str = g[np.where(abs(z-z_str)==delta)]
if isinstance(g_str, list) or isinstance(g_str, np.ndarray):
g_str = g_str[0]
else:
g_str = F*.9999
return g_str
def long_coef(F, Es, z_inflow, z_initial, g_in=1.):
'''
Supplemental function for equation z'(g) after RPS has started.
------
Input:
------
F normalized inflow, xi
Es normalized forced gas outflow, zeta
z_inflow alpha-abundance of inflow, expressed as z'_inflow
z_initial initial alpha-abundance, for z'_initial
g_in initial fraction of gas
------
Returns:
------
Constant term for equation z'(g) after RPS if F!=Es.
'''
return ((g_in/F+(Es-F)/F)*hyp2f1(1,1,F/(F-Es)+1, g_in/(F-Es))+(F-Es)/F+z_inflow-z_initial)
def eq_after_g_abs(F, Es, z_inflow, z_initial, long_coef, g, g_in=1.):
'''
Determines z' = Z/p_eff = (10^[alpha/H])/p_eff for
array of g using equation "After RPS" where F != Es
(where p_eff = p/(1+eta) - effective yield).
------
Input:
F normalized inflow, xi
Es normalized forced gas outflow, zeta
z_inflow alpha-abundance of inflow, expressed as z'_inflow
z_initial initial alpha-abundance, for z'_initial
long_coef constant value (see long_coef() function)
g array of g values
g_in gas at beginning of RPS epoch, g_stripped
------
Returns:
z array of z' values: len(z)=len(g)
'''
return -np.power(g*(g_in-F+Es)/(g_in*(g-F+Es)), F/(Es-F))*long_coef + \
1/F*(g-F+Es)*hyp2f1(1,1,F/(F-Es)+1, g/(F-Es)) + (F-Es)/F+z_inflow
def eq_after_g_equal(F, z_inflow, z_initial, g, g_in):
'''
Determines z' = Z/p_eff = (10^[alpha/H])/p_eff for
array of g using equation "After RPS" where F = Es
(where p_eff = p/(1+eta) - effective yield).
------
Input:
F normalized inflow xi
z_inflow alpha-abundance of the inflow expressed as z'_inflow
z_initial initial alpha-abundance for z'_initial
g array of g values
g_in initial fraction of gas
------
Returns:
z array of z' values: len(z)=len(g).
'''
return np.exp(-F*(1./g-1./g_in))*(z_initial-z_inflow) + \
np.exp(-F/g)*(expi(F/g)-expi(F/g_in)) + z_inflow
def gas_equation(p_eff, F, Es, aH_inflow, aH_initial, aH_str, points = 5001, show=True):
'''
Function to determine g(z').
------
Input:
------
p_eff effective yield in Solar metallicity units
F normalized inflow, xi
Es normalized forced gas outflow, zeta
aH_inflow alpha-abundance of inflow
aH_initial initial alpha-abundance
aH_str alpha-abundance at moment when RPS begins
points length of returned g array
show show the figure of g([alpha/H])
------
Returns:
------
g([alpha/H]) gas as a function of alpha-abundance
'''
# Determine two z'(g) (before and after RPS) separately.
# Then consider only parts which are required to construct final g(z')
arr = 1000
z_inflow = (10**aH_inflow)/p_eff
z_initial = (10**aH_initial)/p_eff
z_str = (10**aH_str)/p_eff
# If xi<1, g will decrease with [alpha/H] for [alpha/H]<[alpha/H]_s
if F<1: loss = True; g_max = 1.
else: loss = False; g_max = 1.9999*F
g = np.linspace(0, g_max, arr)
# Initial gas g_in of equation "After RPS":
g_str = g_rps(F, z_inflow, z_initial, z_str)
# "After RPS" inflow xi:
F2 = 0.01
# Derive long_coef for "After RPS" equation
if F2!=Es:
coef_here = long_coef(F2, Es, z_inflow, z_str, g_str)
# Calculate z'(g) from "Before RPS"
z_before = eq_before_g(F, z_inflow, z_initial, g)
# If xi>1, and if p_eff is very small, g achieves xi at max(z_before) << z_str,
# which from Extreme Inflow model of Larson (1986) balances gas loss with gas inflow
# and makes further increase of alpha-abundance impossible:
if np.max(z_before[np.where(np.isfinite(z_before))])<z_str:
return
z_after = []
# Calculate z'(g) from "After RPS".
for i in range(len(g)):
if F2==Es:
z_after.append(eq_after_g_equal(F2, z_inflow, z_str, g[i], g_str))
else:
z_after.append(eq_after_g_abs(F2, Es, z_inflow, z_str, coef_here, g[i], g_str))
z_before = np.array(z_before)
z_after = np.array(z_after)
z_before, g_before = z_before[np.where(np.isfinite(z_before) & (z_before>0))], g[np.where(np.isfinite(z_before) & (z_before>0))]
aH_before = np.log10(z_before*p_eff)
z_after, g_after = z_after[np.where(np.isfinite(z_after) & (z_after>0)) ], g[np.where(np.isfinite(z_after) & (z_after>0))]
# z'(g) for "After RPS" with Es>F2 is convex function.
# On g(z') plot, let's consider it as two branches, upper and lower.
# for upper branch z' increases with gas decrease
# for lower branch z' decreases with gas decrease
# Conditions for IRPS usage were determined to exclude that (lower) metal-poor branch:
z_max = z_after.max()
z_after0 = z_after
g_after0 = g_after
g_edge = g_after[np.where(z_after==z_max)]
z_after, g_after = z_after[np.where(g_after>g_edge)], g_after[np.where(g_after>g_edge)]
aH_after = np.log10(z_after*p_eff)
# In both cases g increments. If loss==True or Es>F2, [alpha/H] will
# decrease. To recover order arrays must be flipped as follows
if loss:
aH_before, g_before = np.flip(aH_before), np.flip(g_before)
if loss or Es>=F2:
aH_after, g_after = np.flip(aH_after), np.flip(g_after)
g_f = g_before[np.where(aH_before<=aH_str)]
aH_f = aH_before[np.where(aH_before<=aH_str)]
g_f = np.append(g_f, g_after[np.where(aH_after>aH_str)])
aH_f = np.append(aH_f, aH_after[np.where(aH_after>aH_str)])
if show:
plt.plot(aH_before, g_before, color='blue', linestyle='dashed', label='g([alpha/H])\nBefore RPS')
plt.plot(np.log10(z_after0*p_eff), g_after0, color='y', linestyle='dashed', label='After RPS (all)')
plt.plot(aH_after, g_after, color='green', linestyle='solid', lw=0.5, label='After RPS')
plt.ylim([-2.5, F*1.1])
plt.xlim([-5., 0.])
g_f, aH_f = g_f[np.where((np.isfinite(g_f) & np.isfinite(aH_f) & (aH_f<5.)))],\
aH_f[np.where((np.isfinite(g_f) & np.isfinite(aH_f) & (aH_f<5.)))]
aH_min_1 = np.min(aH_before)
aH_max_1 = np.max(aH_after)
# Interpolate for equally spaced [alpha/H] values.
# Absent pieces are fulfilled with g=1, g=0 for beginning and end of
# evolution, if needed
alphah = np.linspace(-5., 0., points)
gas = np.interp(alphah, aH_f, g_f)
gas[np.where(alphah<aH_min_1)] = 1.
if aH_max_1 < 0.:
gas[np.where(alphah>=aH_max_1)] = 0.
if show:
plt.plot(aH_f, g_f, color = 'cyan', linestyle='solid', label = 'Construct' )
plt.plot(np.linspace(-5., 0., points), gas, color = 'k', linestyle='dotted', label = 'Interpolated')
plt.xlabel("[alpha/H]")
plt.ylabel("Gas")
return gas
Part 2: IRPS model
The α-element distribution function (ADF) of IRPS model you can find from a2() function. For the input, the gas points should be predetermined for each of alphaH points from gas_equation function. Keep in mind that alphaH is between -5 and 0 with length points points (argument of gas_equation function)
In [3]:
def a2(p, aH_f, F, E, alphah, gas, aH_str, show=False):
'''
The IRPS model (which sometimes I refer as Accretion-2, A2).
------
Input:
------
p effective yield in Solar metallicity units
aH_f alpha-abundance of inflow
F normalized inflow xi
E normalized forced gas outflow zeta'_s
alphah array of [alpha/H] where the alpha-abundance distribution function should be calculated
gas array of gas points for each [alpha/H]
aH_str alpha-abundance at the moment when RPS begins, aH_s
------
Returns:
------
r([alpha/H]) Model points.
'''
l = None
if isinstance(gas, float):
gas = np.array([gas])
alphah = np.array([alphah])
l = 1.
r = np.zeros(len(gas))
ph1 = np.where(((gas!=0)&(alphah<=aH_str)))[0]
ph2 = np.where(((gas!=0)&(alphah>aH_str)))[0]
r[ph1] = (10**alphah[ph1])*gas[ph1]/p/(1+(10.**aH_f-10.**alphah[ph1])/p*F/gas[ph1])
r[ph2] = (10**alphah[ph2])*gas[ph2]/p/(1+(10.**aH_f-10.**alphah[ph2])/p*0.01/gas[ph2])
r[np.where(r<0.)] = 0.
if isinstance(l, float):
r = r[0]
show = False
if show:
p = plt.plot(alphah, gas, ls = 'dotted', label='Gas')
c = p[0].get_color()
plt.plot(alphah, r, c=c, ls = 'solid', lw=0.75, label = "ADF: \u03BE"+f'={F}, '+'[alpha/H]$_s$'+f'={aH_str}, '+'\u03B6'+f'={E}')
plt.xlabel("[alpha/H]")
plt.ylabel("ADF")
return r
Examples: fitting [α/H]str, peff, ζ/ξ for IRPS ADF
Prompt values for g([$\alpha$/H]), $\alpha$-elements Distribution Function (ADF)
In [4]:
alphah = np.linspace(-5., 0., 5001)
p = 0.0222; F = 30.2; E = 1.4; aH_str = -1.73
format_ticks()
gas2 = gas_equation(0.05, 30.01, 15., -5, -5, aH_str=-1.85, points=5001, show=True);
plt.legend(bbox_to_anchor=(1.04, 1), loc="upper left");
plt.title(r"$\mathbf{Figure~1.}$ Inside of the $\mathit{gas \_ equation}$ function: gas-abundance relation.");
plt.show()
format_ticks()
a2(p, -5., F, E, alphah, gas_equation(p, F, E, -5, -5, aH_str, 5001, show=False), aH_str = aH_str, show=True);
plt.legend(bbox_to_anchor=(1.04, 1), loc="upper left");
plt.title(r"$\mathbf{Figure~2.}$ Inside of the $\mathit{a2}$ function: gas-abundance relation and ADF.");
In [5]:
# Try it
Example #1: $\alpha$-abundance of ram pressure stripping [$\alpha$/H]$_s$
In [6]:
p = 0.75; F = 10; E = 10.5; aH_str = -1.3; alphah = np.linspace(-5., 0., 5001)
format_ticks()
for fs in range(-20, -11, 2):
a2(p, -5., F, E, alphah, gas_equation(p, F, E, -5, -5, fs/10., show=False), aH_str = fs/10.,\
show=True)
plt.legend(bbox_to_anchor=(1.04, 1), loc="upper left")
plt.title(r"$\mathbf{Figure~3.}$ IRPS ADFs for different [$\alpha$/H]$_{str}$ with other parameters fixed");
plt.show()
In [7]:
# Try it
Example #2: inflow $\xi$
In [8]:
p = 0.75; F = 10; E = 10.5; aH_str = -1.3; alphah = np.linspace(-5., 0., 5001)
format_ticks()
for f in range(11, 300, 60):
a2(p, -5., f/10., E, alphah, gas_equation(p, f/10., E, -5, -5, aH_str, show=False), aH_str = aH_str,\
show=True)
plt.legend(bbox_to_anchor=(1.04, 1), loc="upper left")
plt.title(r"$\mathbf{Figure~4.}$ IRPS ADFs for different inflow $\xi$ ($\mathit{system's~gas~inflow}$) with other parameters fixed");
plt.show()
In [9]:
# Try it too
Example #3.1: terminal wind rate $\zeta$
In [10]:
p = 0.75; F = 10; E = 10.5; aH_str = -1.3; alphah = np.linspace(-5., 0., 5001)
format_ticks()
for e in range(11, 300, 60):
a2(p, -5., F, e/10., alphah, gas_equation(p, F, e/10., -5, -5, aH_str, show=False), aH_str = aH_str,\
show=True)
plt.legend(bbox_to_anchor=(1.04, 1), loc="upper left")
plt.title(r"$\mathbf{Figure~5.}$ IRPS ADFs for different $\zeta$ ($\mathit{forced~gas~loss}$) with other parameters fixed");
plt.show()
In [11]:
# Try it
Example #3.2: terminal wind rate $\zeta$, the second version
In [12]:
# Not required, but you can try it to see the shape of ADF (alpha/H, zeta)
plotly.offline.init_notebook_mode()
p = 0.75; F = 100.; aH_str = -1.3;
alphah = np.linspace(-5., 0., 50)
Es = np.linspace(2., 190, 50)
mdf = np.zeros((len(alphah), len(Es)))
for j in range(len(Es)):
gas_loc = gas_equation(p, F, Es[j], -5, -5, aH_str, len(alphah), False)
for i in range(len(alphah)):
mdf[i][j] = a2(p, -5., F, Es[j], alphah[i], gas_loc[i], aH_str, False)
X = Es
Y = alphah
Z = np.reshape(mdf, (len(alphah), len(Es)))
layout = go.Layout(margin={'l': 0, 'r': 0, 'b': 0, 't': 0})
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z = dict(show=True, usecolormap=True, highlightcolor="aliceblue", project_z=True),
contours_x = dict(show=True, usecolormap=True, highlightcolor="aliceblue", project_x=True))
fig.update_layout(title=r'Figure 6. Similar to Figure 5, but in 3D.',
autosize=False,
scene_camera_eye=dict(x=1.87, y=1.2, z=.8),
scene = dict(
xaxis_title='zeta',
yaxis_title='[alpha/H]',
zaxis_title='ADF',
zaxis = dict(nticks=15, range=[-1,mdf.max()],),),
width=650, height=650,
margin=dict(l=65, r=50, b=65, t=90))
fig.show()
About this notebook
In case of question, contact kkvasova@nd.edu
Citations
If you use numpy, scipy, matplotlib, plotly for published research, please cite the authors. Follow these links for more information about citing:
In [ ]:
In [ ]: